Intro

Fit density (also referred to as CPUE) model with environmental predictors and use that to calculate weighted mean dissolved oxygen, temperature and depth of Baltic cod

# Load libraries, install if needed
library(tidyverse); theme_set(theme_light(base_size = 11))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(qwraps2) 
library(wesanderson)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)

# To load entire cache in interactive r session, do: 
# qwraps2::lazyload_cache_dir(path = "R/analysis/cpue_model_cache/html")

For maps

# Specify map ranges
ymin = 52; ymax = 60; xmin = 10; xmax = 24

map_data <- rnaturalearth::ne_countries(
  scale = "medium",
  returnclass = "sf", continent = "europe")

# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
  st_crop(map_data,
          c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))

# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)

# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
      axis.text.x = element_text(angle = 90),
      axis.text = element_text(size = 8),
      legend.text = element_text(size = 8),
      legend.title = element_text(size = 8),
      legend.position = "bottom",
      legend.key.height = unit(0.2, "cm"),
      legend.margin = margin(0, 0, 0, 0),
      legend.box.margin = margin(-5, -5, -5, -5),
      strip.text = element_text(size = 8, colour = 'gray10', margin = margin()),
      strip.background = element_rect(fill = "grey95")
      )
}

# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
  theme_light(base_size = 10, base_family = "") +
    theme(
        axis.text.x = element_text(angle = 90),
        axis.text = element_text(size = 6),
        strip.text = element_text(size = 8, colour = 'gray10', margin = margin()),
        strip.background = element_rect(fill = "gray95"),
        legend.position = c(0.7, 0.02),
        legend.direction = "horizontal"
      )
}

# Make default base map plot
xmin2 <- 303000; xmax2 <- 916000; xrange <- xmax2 - xmin2
ymin2 <- 5980000; ymax2 <- 6450000; yrange <- ymax2 - ymin2

plot_map_raster <- 
ggplot(swe_coast_proj) + 
  xlim(xmin2, xmax2) +
  ylim(ymin2, ymax2) +
  labs(x = "Longitude", y = "Latitude") +
  geom_sf(size = 0.3) + 
  theme_plot()

plot_map_raster_labels <- 
  plot_map_raster + 
  annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = 1.9, angle = 75) +
  annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = 1.9) +
  annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = 1.9, angle = 75) +
  annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = 1.9, angle = 75)

Read data

d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cpue_q_1_4.csv")

# Calculate standardized variables
d <- d %>% 
  filter(year < 2020 & year > 1992 & quarter == 4) %>% 
  mutate(oxy_sc = oxy,
         temp_sc = temp,
         depth_sc = depth,
         ) %>%
  mutate_at(c("oxy_sc", "temp_sc", "depth_sc"),
            ~(scale(.) %>% as.vector)) %>% 
  mutate(year = as.integer(year)) %>% 
  drop_na(oxy, depth, temp)

plot_map_raster +
  geom_point(data = d, aes(X*1000, Y*1000))

Read the prediction grids

pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
pred_grid <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")

pred_grid <- bind_rows(pred_grid1, pred_grid)

# Standardize data with respect to data grid:
pred_grid <- pred_grid %>%
  mutate(year = as.integer(year)) %>% 
  filter(year %in% c(unique(d$year))) %>% 
  mutate(depth_sc = (depth - mean(d$depth))/sd(d$depth),
         temp_sc = (temp - mean(d$temp))/sd(d$temp),
         oxy_sc = (oxy - mean(d$oxy))/sd(d$oxy)) %>% # Need to scale these to the mean and sd in the data!
  drop_na(oxy, depth, temp)

Make spde mesh

spde <- make_mesh(d, xy_cols = c("X", "Y"),
                  n_knots = 200, 
                  type = "kmeans", seed = 42)

# Plot and save spde
png(file = "figures/supp/density/spde.png", units = "in", width = 6.5, height = 6.5, res = 300)
plot(spde)
dev.off()

Fit the density model

# Depth spline + oxy spline
# Takes about 30 minutes
m <- sdmTMB(density ~ 0 + as.factor(year) + s(depth_sc) + s(oxy_sc) + s(temp_sc),
            data = d,
            mesh = spde,
            family = tweedie(link = "log"),
            spatiotemporal = "AR1",
            spatial = "on",
            time = "year",
            reml = FALSE,
            control = sdmTMBcontrol(newton_steps = 1))
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0163339061734842.

tidy(m, conf.int = TRUE)

Sanity check

sanity(m)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> x `thetaf` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_phi` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large

m_1 <- run_extra_optimization(m, nlminb_loops = 0, newton_loops = 1)

sanity(m_1)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> ✓ No gradients with respect to fixed effects are >= 0.001
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large

Check QQ plot and residuals

mcmc_res <- residuals(m_1, type = "mle-mcmc", mcmc_iter = 5000, mcmc_warmup = 2000,
                      stan_args = list(control = list(max_treedepth = 11)))
#> Warning: There were 5 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 11. See
#> http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is 1.06, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#bulk-ess

png(file = "figures/supp/density/qq.png", units = "in", width = 6.5, height = 6.5, res = 300)
qqnorm(mcmc_res);qqline(mcmc_res)
dev.off()

d$residuals <- mcmc_res[, 1]
# Residuals on map
plot_map_raster +
  geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = residuals)) +
  scale_colour_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  theme(strip.text = element_text(size = 8, colour = 'black', margin = margin()),
        strip.background = element_rect(fill = "grey90"))
#> Warning: Removed 8 rows containing missing values (geom_point).


ggsave("figures/supp/density/spatial_resid.png", width = 6.5, height = 6.5, dpi = 600)
#> Warning: Removed 8 rows containing missing values (geom_point).

Check the AR1 parameter (rho is ar_phi on the -1 to 1 scale):

tidy(m_1, effects = "ran_pars", conf.int = TRUE)

Predict on grid

predict_mcod <- predict(m_1, newdata = pred_grid)

# Save prediction as csv (smaller than the model object...)
write.csv(predict_mcod, "output/predict_mcod_density.csv")

Plot predictions on map

# Plot predicted density and random effects
plot_map_raster +
  geom_raster(data = filter(predict_mcod, depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
  scale_fill_viridis_c(trans = "sqrt") +
  facet_wrap(~ year, ncol = 5) +
  labs(fill = expression(kg/km^2)) +
  ggtitle("Predicted density (fixed + random)")
#> filter: removed 36,936 rows (15%), 212,517 rows remaining

 
# I save this below instead to better match the main figure 
#ggsave("figures/supp/density/est_map.png", width = 6.5, height = 6.5, dpi = 600)

# Plot spatiotemporal random effect
plot_map_raster +
  geom_raster(data = predict_mcod, aes(x = X * 1000, y = Y * 1000, fill = epsilon_st)) +
  scale_fill_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  ggtitle("Spatiotemporal random effects")


ggsave("figures/supp/density/epsilon_st_map.png", width = 6.5, height = 6.5, dpi = 600)

# Plot spatial random effect
plot_map_raster +
  geom_raster(data = filter(predict_mcod, year == 1999), aes(x = X * 1000, y = Y * 1000, fill = omega_s)) +
  scale_fill_gradient2() +
  facet_wrap(~ year, ncol = 5) +
  theme_plot() +
  ggtitle("Spatial random effects")
#> filter: removed 240,214 rows (96%), 9,239 rows remaining

  
ggsave("figures/supp/density/omega_s_map.png", width = 6.5, height = 6.5, dpi = 600)

Annual indices using get_index_sims method to avoid slow bias correction (not done above!)

# preds_cod_sim <- predict(m_1, newdata = pred_grid, sims = 100)
# x <- get_index_sims(preds_cod_ave_sim, area = rep(2*2, nrow(preds_cod_ave_sim)))
# 
# ggplot(x, aes(year, est/1000, ymin = lwr/1000, ymax = upr/1000)) +
#     geom_line() +
#     geom_ribbon(alpha = 0.4)
# x_sims <- get_index_sims(preds_cod_ave_sim, return_sims = TRUE)
# ggplot(x_sims, aes(as.factor(year), log(.value))) +
#     geom_violin()

###
# Now get the total index by specifying the area of a grid cell
# Total prediction 
preds_cod_sim <- predict(m_1, newdata = pred_grid, sims = 100)
#> Warning: The `sims` argument of `predict.sdmTMB()` is deprecated as of sdmTMB 0.0.21.
#> Please use the `nsim` argument instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

# Total prediction (omitting SD 24!)
preds_cod25_28_sim <- predict(m_1, newdata = filter(pred_grid, !sub_div == 24), sims = 100)
#> filter: removed 38,421 rows (15%), 211,032 rows remaining

# SD 24-25
preds_cod24_25_sim <- predict(m_1, newdata = filter(pred_grid, sub_div %in% c(24, 25)), sims = 100)
#> filter: removed 140,562 rows (56%), 108,891 rows remaining

# SD 24
preds_cod24_sim <- predict(m_1, newdata = filter(pred_grid, sub_div == 24), sims = 100)
#> filter: removed 211,032 rows (85%), 38,421 rows remaining

# SD 25
preds_cod25_sim <- predict(m_1, newdata = filter(pred_grid, sub_div == 25), sims = 100)
#> filter: removed 178,983 rows (72%), 70,470 rows remaining

# SD 26
preds_cod26_sim <- predict(m_1, newdata = filter(pred_grid, sub_div == 26), sims = 100)
#> filter: removed 184,410 rows (74%), 65,043 rows remaining

# SD 27
preds_cod27_sim <- predict(m_1, newdata = filter(pred_grid, sub_div == 27), sims = 100)
#> filter: removed 224,505 rows (90%), 24,948 rows remaining

# SD 28
preds_cod28_sim <- predict(m_1, newdata = filter(pred_grid, sub_div == 28), sims = 100)
#> filter: removed 198,882 rows (80%), 50,571 rows remaining

# Now calculate the CPUE index (average)
index_sim <- get_index_sims(preds_cod_sim, area = rep(2*2, nrow(preds_cod_sim)))
index24_25_sim <- get_index_sims(preds_cod24_25_sim, area = rep(2*2, nrow(preds_cod24_25_sim)))
index25_28_sim <- get_index_sims(preds_cod25_28_sim, area = rep(2*2, nrow(preds_cod25_28_sim)))
index24_sim <- get_index_sims(preds_cod24_sim, area = rep(2*2, nrow(preds_cod24_sim)))
index25_sim <- get_index_sims(preds_cod25_sim, area = rep(2*2, nrow(preds_cod25_sim)))
index26_sim <- get_index_sims(preds_cod26_sim, area = rep(2*2, nrow(preds_cod26_sim)))
index27_sim <- get_index_sims(preds_cod27_sim, area = rep(2*2, nrow(preds_cod27_sim)))
index28_sim <- get_index_sims(preds_cod28_sim, area = rep(2*2, nrow(preds_cod28_sim)))

# Join indices
index_sim <- index_sim %>% mutate(sub_div = "Total")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
index24_sim <- index24_sim %>% mutate(sub_div = "24")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
index25_sim <- index25_sim %>% mutate(sub_div = "25")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
index26_sim <- index26_sim %>% mutate(sub_div = "26")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
index27_sim <- index27_sim %>% mutate(sub_div = "27")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
index28_sim <- index28_sim %>% mutate(sub_div = "28")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA

div_index_sim <- bind_rows(index_sim, index24_sim, index25_sim, index26_sim, index27_sim, index28_sim)

Compare the above biomass-index with one where I filter deep depths

# Now do the same but excluding the areas not sampled in the pred grid
ncells_sub <- filter(pred_grid, year == max(pred_grid$year) & depth > 20 & depth < 120) %>% nrow()
#> filter: removed 242,779 rows (97%), 6,674 rows remaining
pred_grid_sub <- filter(pred_grid, depth > 20 & depth < 120)
#> filter: removed 69,255 rows (28%), 180,198 rows remaining

pred_avg_sim_sub <- predict(m_1, newdata = pred_grid_sub, nsim = 100)

index_avg_sim_sub <- get_index_sims(pred_avg_sim_sub, area = rep(2*2, nrow(pred_avg_sim_sub)))

# Combine and plot & compare
index_avg_sim_comp <- bind_rows(mutate(index_sim, area = "full"),
                                mutate(index_avg_sim_sub, area = "subset")) %>% 
  mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> mutate: new variable 'est_t' (double) with 54 unique values and 0% NA
#>         new variable 'lwr_t' (double) with 54 unique values and 0% NA
#>         new variable 'upr_t' (double) with 54 unique values and 0% NA

ggplot(index_avg_sim_comp, aes(year, est_t, ymin = lwr_t, ymax = upr_t, color = area)) +
  geom_line() +
  geom_ribbon(alpha = 0.4)

Combined plot with biomass index by sd, total and predictions on map for two years

div_index_sim <- div_index_sim %>% mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000) 

ind_plot_sd <- 
  ggplot() +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  geom_line(data = filter(div_index_sim, !sub_div == "Total"), aes(year, est_t, color = sub_div)) + 
  geom_ribbon(data = filter(div_index_sim, !sub_div == "Total"), aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "ICES Sub division") +
  guides(color = "none") + 
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme_plot() + 
  theme(axis.text.x = element_text(angle = 0),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL

ind_plot_tot <- 
  ggplot() +
  geom_line(data = filter(div_index_sim, sub_div == "Total"),
            aes(year, est_t)) +
  geom_ribbon(data = filter(div_index_sim, sub_div == "Total"),
              aes(year, ymin = lwr_t, ymax = upr_t, fill = "Total"), alpha = 0.4) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_fill_manual(values = "grey40") + 
  labs(x = "Year", y = "Tonnes", fill = "") +
  theme_plot() + 
  theme(axis.text.x = element_text(angle = 0),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.2, 0.15, 0), "cm")) +
  NULL

cpue_map_94_18 <- plot_map_raster_labels +
  geom_raster(data = filter(predict_mcod, year %in% c("1994", "2018") & depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
  scale_fill_viridis_c(trans = "sqrt") +
  facet_wrap(~ year, ncol = 2) +
  labs(fill = expression(kg/km^2)) +
  theme_plot() +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
        legend.text = element_text(size = 5, angle = 90),
        legend.position = "bottom") +
  NULL

(ind_plot_sd | ind_plot_tot) / cpue_map_94_18 + plot_annotation(tag_levels = 'A')


ggsave("figures/density.png", width = 6.5, height = 6.5, dpi = 600)

# All years cpue
plot_map_raster +
  geom_raster(data = filter(predict_mcod, depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
  scale_fill_viridis_c(trans = "sqrt") +
  facet_wrap(~ year, ncol = 5) +
  labs(fill = expression(kg/km^2)) +
  theme_plot() +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
        legend.text = element_text(size = 5, angle = 90),
        legend.position = "bottom") +
  NULL


ggsave("figures/supp/density/all_years_condition_map_covars.png", width = 6.5, height = 6.5, dpi = 600)

Biomass weighted quantiles of depth and oxygen

# Plot bathymetry
bathy_plot <- plot_map_raster_labels +
  geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = depth)) +
  scale_fill_viridis(option = "A", direction = -1) +
  labs(fill = "Depth") +
  theme_plot() +
  theme(legend.position = "bottom",
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  NULL
#> filter: removed 240,214 rows (96%), 9,239 rows remaining
  
# Now calculate quantiles of the biomass-weighted values
pal <- c(rev(brewer.pal(n = 8, name = "Dark2")), "gray20")

wm_depth <- predict_mcod %>%
  group_by(year) %>%
  summarise("Median" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.5)),
            "1st quartile" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.25)),
            "3rd quartile" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.75))) %>% 
  pivot_longer(cols = c("Median", "1st quartile", "3rd quartile"),
               names_to = "series", values_to = "depth") %>% 
  group_by(series) %>%  # standardize within for easy plotting
  mutate(depth_ct = depth - mean(depth)) %>% 
  ungroup()
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> pivot_longer: reorganized (Median, 1st quartile, 3rd quartile) into (series, depth) [was 27x4, now 81x3]
#> group_by: one grouping variable (series)
#> mutate (grouped): new variable 'depth_ct' (double) with 72 unique values and 0% NA
#> ungroup: no grouping variables
  
plot_w_dep <- ggplot(wm_depth, aes(year, depth, color = series, group = series, fill = series)) +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 1) +
  geom_point(size = 1.5, alpha = 0.8) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_manual(values = pal, name = "Quantiles") +
  scale_fill_manual(values = pal) +
  guides(fill = "none", color = guide_legend(nrow = 2, title.position = "top", title.hjust = 0.5)) +
  labs(y = "Depth [m]", x = "Year", color = "") +
  scale_y_reverse() + 
  theme_plot() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 0),
        legend.position = "bottom",
        legend.key.size = unit(0.6, "cm"),
        legend.text = element_text(size = 8),
        plot.margin = unit(c(0, 0.4, 0, 0), "cm")) +
  NULL

# Plot oxygen in 2006
oxy_plot <- plot_map_raster_labels +
  geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = oxy)) +
  scale_fill_viridis() +
  labs(fill = expression(paste("O" [2], " [ml/L]", sep = ""))) +
  theme_plot() +
  theme(legend.position = "bottom",
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  NULL
#> filter: removed 240,214 rows (96%), 9,239 rows remaining

# Plot biomass-weighted oxygen
# Weighted total
wm_oxy <- predict_mcod %>%
  group_by(year) %>%
  summarise("Median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)),
            "1st quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.25)),
            "3rd quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.75))) %>% 
  pivot_longer(cols = c("Median", "1st quartile", "3rd quartile"),
               names_to = "series", values_to = "oxy") %>% 
  ungroup() %>% 
  group_by(series) %>% # standardize within for easy plotting
  mutate(oxy_ct = oxy - mean(oxy)) %>%
  ungroup()
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> pivot_longer: reorganized (Median, 1st quartile, 3rd quartile) into (series, oxy) [was 27x4, now 81x3]
#> ungroup: no grouping variables
#> group_by: one grouping variable (series)
#> mutate (grouped): new variable 'oxy_ct' (double) with 81 unique values and 0% NA
#> ungroup: no grouping variables


# Find which depths to calculate average oxygen on

wm_depth %>% filter(series == "1st quartile") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped
wm_depth %>% filter(series == "3rd quartile") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped

annual_oxy <- pred_grid %>%
  drop_na(oxy) %>%
  filter(depth > 29 & depth < 61) %>% ###
  #filter(depth > 40 & depth < 60) %>% ###
  group_by(year) %>%
  summarize(median_oxy = median(oxy)) %>%
  ungroup() %>% 
  rename("oxy" = "median_oxy") %>% 
  mutate(series = "Median (env.)",
         oxy_ct = oxy - mean(oxy))
#> drop_na: no rows removed
#> filter: removed 179,847 rows (72%), 69,606 rows remaining
#> group_by: one grouping variable (year)
#> summarize: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> rename: renamed one variable (oxy)
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#>         new variable 'oxy_ct' (double) with 27 unique values and 0% NA
  
oxygen_series <- bind_rows(wm_oxy, annual_oxy)

plot_w_oxy <-
ggplot(oxygen_series, aes(year, oxy, color = series, group = series, fill = series, linetype = series)) +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 1) +
  geom_point(size = 1.5, alpha = 0.8) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_linetype_manual(values = c(1,1,1,2)) +
  scale_color_manual(values = pal, name = "Quantiles") +
  scale_fill_manual(values = pal) +
  guides(fill = "none", linetype = "none", color = guide_legend(nrow = 2, title.position = "top", title.hjust = 0.5)) +
  labs(y = expression(paste("O" [2], " [ml/L]", sep = "")), x = "Year",
       color = "") +
  theme_plot() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 0),
        legend.position = "bottom",
        legend.key.size = unit(0.6, "cm"),
        legend.text = element_text(size = 8),
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  NULL

plot_w_oxy


(bathy_plot | oxy_plot) / (plot_w_dep | plot_w_oxy) +
  plot_annotation(tag_levels = 'A') 


ggsave("figures/weighted_quantiles.png", width = 6.5, height = 6.5, dpi = 600)

# Now see how the weighted mean oxygen differs from the average oxygen at the range of depths used by cod
wm_depth %>% group_by(series) %>% summarise(mean_depth = mean(depth))
#> group_by: one grouping variable (series)
#> summarise: now 3 rows and 2 columns, ungrouped

wm_oxy2 <- wm_oxy %>% filter(series == "Median") %>% mutate(series = "Biomass-weighted")
#> filter: removed 54 rows (67%), 27 rows remaining
#> mutate: changed 27 values (100%) of 'series' (0 new NA)

pred_grid_ave_oxy <- pred_grid %>%
  filter(depth < 60 & depth > 30) %>% 
  group_by(year) %>% 
  summarise(oxy = mean(oxy)) %>% 
  mutate(series = "Mean oxygen at depth 30-60 m")
#> filter: removed 185,247 rows (74%), 64,206 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
#> mutate: new variable 'series' (character) with one unique value and 0% NA

combi <- bind_rows(wm_oxy2, pred_grid_ave_oxy)

ggplot(combi, aes(year, oxy, color = series)) +
  geom_point() +
  labs(y = expression(paste("O" [2], " [ml/L]", sep = "")), x = "Year",
       color = "") +
  scale_color_manual(values = pal[c(1,3)], name = "") +
  stat_smooth(method = "lm") +
  theme_plot()
#> `geom_smooth()` using formula 'y ~ x'


ggsave("figures/supp/density/weighted_oxy_vs_depth_oxy.png", width = 6.5, height = 6.5, dpi = 600)
#> `geom_smooth()` using formula 'y ~ x'

Weighted oxygen and depth by year and sub div

# Oxygen
wm_oxy_sd <- predict_mcod %>%
  group_by(year, sub_div) %>%
  summarise("weighted_median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)))
#> group_by: 2 grouping variables (year, sub_div)
#> summarise: now 135 rows and 3 columns, one group variable remaining (year)

ggplot(wm_oxy_sd, aes(year, weighted_median, color = factor(sub_div))) +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 4), size = 1) +
  geom_point(size = 2.5, alpha = 0.8) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = "none", fill = "none", linetype = "none") +
  labs(y = expression(paste("O" [2], " [ml/L]", sep = "")), x = "Year",
       color = "") +
  theme_plot() +
  facet_wrap(~sub_div) +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 0),
        legend.position = "bottom",
        legend.key.size = unit(0.6, "cm"),
        legend.text = element_text(size = 8),
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  NULL


ggsave("figures/supp/density/weighted_oxygen_subdiv.png", width = 6.5, height = 6.5, dpi = 600)

# # Calculate weighted mean and quantiles for sd 25 only 
predict_mcod %>%
  filter(sub_div == 25 & year >= 2015) %>%
  #group_by(year) %>%
  summarise("Median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)),
            "1st quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.25)),
            "3rd quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.75))) %>%
  as.data.frame()
#> filter: removed 236,403 rows (95%), 13,050 rows remaining
#> summarise: now one row and 3 columns, ungrouped

# Now compare delta change in oxygen with delta change in condition
cond_delta_sd <- read.csv("output/div_index_delta_cond.csv") %>%
  mutate(sub_div = as.character(sub_div)) %>%
  dplyr::select(-X) %>%
  filter(!sub_div == "Total")
#> mutate: no changes
#> filter: removed one row (17%), 5 rows remaining

wm_oxy_sd_delta <- wm_oxy_sd %>%
  group_by(sub_div) %>%
  filter(year %in% c(1993, 1994, 1995, 2017, 2018, 2019)) %>%
  dplyr::select(year, weighted_median, sub_div) %>%
  pivot_wider(1:3, names_from = year, values_from = weighted_median) %>%
  mutate(start_oxy = mean(`1993`, `1994`, `1995`),
         end_oxy = mean(`2017`, `2018`, `2019`),
         delta_oxy = end_oxy - start_oxy) %>%
  ungroup() %>%
  dplyr::select(sub_div, delta_oxy) %>%
  mutate(sub_div = as.character(sub_div))
#> group_by: one grouping variable (sub_div)
#> filter (grouped): removed 105 rows (78%), 30 rows remaining
#> pivot_wider: reorganized (year, weighted_median) into (1993, 1994, 1995, 2017, 2018, …) [was 30x3, now 5x7]
#> mutate (grouped): new variable 'start_oxy' (double) with 5 unique values and 0% NA
#>                   new variable 'end_oxy' (double) with 5 unique values and 0% NA
#>                   new variable 'delta_oxy' (double) with 5 unique values and 0% NA
#> ungroup: no grouping variables
#> mutate: converted 'sub_div' from double to character (0 new NA)


# delta_dat <- cond_delta_sd %>% left_join(wm_oxy_sd_delta)
# 
# ggplot(delta_dat, aes(delta_oxy, delta_cond, color = sub_div)) +
#   geom_point(size = 3) +
#   geom_vline(xintercept = 0, color = "red", linetype = 2) +
#   geom_hline(yintercept = 0, color = "red", linetype = 2)
#   NULL


# Depth
wm_dep_sd <- predict_mcod %>%
  group_by(year, sub_div) %>%
  summarise("weighted_median" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.5)))
#> group_by: 2 grouping variables (year, sub_div)
#> summarise: now 135 rows and 3 columns, one group variable remaining (year)

ggplot(wm_dep_sd, aes(year, weighted_median, color = factor(sub_div))) +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 4), size = 1) +
  geom_point(size = 2.5, alpha = 0.8) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = "none", fill = "none", linetype = "none") +
  labs(y = "Depth [m]", x = "Year", color = "") +
  theme_plot() +
  facet_wrap(~sub_div) +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 0),
        legend.position = "bottom",
        legend.key.size = unit(0.6, "cm"),
        legend.text = element_text(size = 8),
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  NULL


ggsave("figures/supp/density/weighted_depth_subdiv.png", width = 6.5, height = 6.5, dpi = 600)

Plot biomass-weighted temperature

# Plot temperature in 2006
temp_plot <- plot_map_raster +
  geom_raster(data = filter(predict_mcod, year == "1999"), aes(x = X * 1000, y = Y * 1000, fill = temp)) +
  scale_fill_viridis(option = "inferno", direction = -1) +
  labs(fill = "Temperature [°C]") +
  theme_plot() +
  theme(legend.position = "bottom",
        plot.margin = unit(c(0, 0.1, 0, 0.2), "cm")) +
  NULL

annual_temp <- pred_grid %>%
  drop_na(temp) %>%
  group_by(year) %>%
  summarize(median_temp = median(temp)) %>%
  ungroup() %>%
  rename("temp" = "median_temp") %>%
  mutate(series = "Median (env.)",
         temp_ct = temp - mean(temp))

# Weighted total
wm_temp <- predict_mcod %>%
  group_by(year) %>%
  summarise("Median" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.5)),
            "1st quartile" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.1)),
            "3rd quartile" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.9))) %>%
  pivot_longer(cols = c("Median", "1st quartile", "3rd quartile"),
               names_to = "series", values_to = "temp") %>%
  ungroup() %>%
  group_by(series) %>% # standardize within for easy plotting
  mutate(temp_ct = temp - mean(temp)) %>%
  ungroup()

temp_series <- bind_rows(wm_temp, annual_temp)

plot_w_temp <- ggplot(temp_series, aes(year, temp, color = series, group = series, fill = series, linetype = series)) +
  stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 1) +
  geom_point(size = 2.5, alpha = 0.8, color = "white", shape = 21) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_linetype_manual(values = c(1,1,1,2)) +
  scale_color_manual(values = pal, name = "Quantiles") +
  scale_fill_manual(values = pal) +
  guides(fill = "none", linetype = "none", color = guide_legend(nrow = 2, title.position = "top", title.hjust = 0.5)) +
  labs(y = "Temperature [°C]", x = "Year",
       color = "") +
  theme_plot() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 0),
        legend.position = "bottom",
        legend.key.size = unit(0.6, "cm"),
        legend.text = element_text(size = 7),
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  NULL

temp_plot + plot_w_temp +
  plot_annotation(tag_levels = 'A') + plot_layout(ncol = 2)


ggsave("figures/supp/density/weighted_temp.png", width = 6.5, height = 8.5, dpi = 600)
annual_temp2 <- pred_grid %>%
  drop_na(temp) %>%
  group_by(year) %>%
  summarize(value = mean(temp)) %>%
  mutate(series = "Mean (env.)",
         variable = "Temperature",
         species = "Environment")
#> drop_na: no rows removed
#> group_by: one grouping variable (year)
#> summarize: now 27 rows and 2 columns, ungrouped
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#>         new variable 'variable' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA

annual_oxy2 <- pred_grid %>%
  drop_na(oxy) %>%
  group_by(year) %>%
  summarize(value = mean(oxy)) %>%
  ungroup() %>%
  mutate(series = "Mean (env.)",
         variable = "Oxygen",
         species = "Environment")
#> drop_na: no rows removed
#> group_by: one grouping variable (year)
#> summarize: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#>         new variable 'variable' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA

# Weighted total
wm_temp2_cod <- predict_mcod %>%
  filter(exp(est) < 3000) %>% 
  group_by(year) %>%
  summarise(value = weighted.mean(x = temp, w = exp(est))) %>%
  ungroup() %>% 
  mutate(series = "Mean (weighted)",
         variable = "Temperature",
         species = "Cod")
#> filter: removed 714 rows (<1%), 248,739 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#>         new variable 'variable' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA

wm_temp2_fle <- predict_mcod %>%
  filter(density_fle < 2000) %>% 
  group_by(year) %>%
  summarise(value = weighted.mean(x = temp, w = density_fle)) %>%
  ungroup() %>% 
  mutate(series = "Mean (weighted)",
         variable = "Temperature",
         species = "Flounder")
#> filter: removed 1,291 rows (1%), 248,162 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#>         new variable 'variable' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA

wm_temp2_spr <- predict_mcod %>%
  group_by(year) %>%
  drop_na(biomass_spr) %>% 
  summarise(value = weighted.mean(x = temp, w = biomass_spr)) %>%
  ungroup() %>% 
  mutate(series = "Mean (weighted)",
         variable = "Temperature",
         species = "Sprat")
#> group_by: one grouping variable (year)
#> drop_na (grouped): removed 21,502 rows (9%), 227,951 rows remaining
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#>         new variable 'variable' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA

wm_temp2_her <- predict_mcod %>%
  group_by(year) %>%
  drop_na(biomass_her) %>% 
  summarise(value = weighted.mean(x = temp, w = biomass_her)) %>%
  ungroup() %>% 
  mutate(series = "Mean (weighted)",
         variable = "Temperature",
         species = "Herring")
#> group_by: one grouping variable (year)
#> drop_na (grouped): removed 21,502 rows (9%), 227,951 rows remaining
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#>         new variable 'variable' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA

wm_oxy2_cod <- predict_mcod %>%
  filter(exp(est) < 3000) %>% 
  group_by(year) %>%
  summarise(value = weighted.mean(x = oxy, w = exp(est))) %>%
  ungroup() %>% 
  mutate(series = "Mean (weighted)",
         variable = "Oxygen",
         species = "Cod")
#> filter: removed 714 rows (<1%), 248,739 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#>         new variable 'variable' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA

wm_oxy2_fle <- predict_mcod %>%
  filter(density_fle < 2000) %>% 
  group_by(year) %>%
  summarise(value = weighted.mean(x = oxy, w = density_fle)) %>%
  ungroup() %>% 
  mutate(series = "Mean (weighted)",
         variable = "Oxygen",
         species = "Flounder")
#> filter: removed 1,291 rows (1%), 248,162 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#>         new variable 'variable' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA

mean_oxy_temp <- bind_rows(annual_temp2, annual_oxy2, wm_temp2_cod, wm_temp2_fle, wm_oxy2_cod, wm_oxy2_fle, wm_temp2_her, wm_temp2_spr)


# coords
wm_coord_cod <- predict_mcod %>%
  filter(exp(est) < 3000) %>% 
  group_by(year) %>%
  summarise(wlat = weighted.mean(x = lat, w = exp(est)),
            wlon = weighted.mean(x = lon, w = exp(est))) %>%
  ungroup() %>% 
  mutate(species = "Cod")
#> filter: removed 714 rows (<1%), 248,739 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'species' (character) with one unique value and 0% NA

wm_coord_fle <- predict_mcod %>%
  filter(density_fle < 2000) %>% 
  group_by(year) %>%
  summarise(wlat = weighted.mean(x = lat, w = density_fle),
            wlon = weighted.mean(x = lon, w = density_fle)) %>%
  ungroup() %>% 
  mutate(species = "Flounder")
#> filter: removed 1,291 rows (1%), 248,162 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'species' (character) with one unique value and 0% NA

wm_coord_spr <- predict_mcod %>%
  group_by(year) %>%
  drop_na(biomass_spr) %>% 
  summarise(wlat = weighted.mean(x = lat, w = biomass_spr),
            wlon = weighted.mean(x = lon, w = biomass_spr)) %>%
  ungroup() %>% 
  mutate(species = "Sprat")
#> group_by: one grouping variable (year)
#> drop_na (grouped): removed 21,502 rows (9%), 227,951 rows remaining
#> summarise: now 27 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'species' (character) with one unique value and 0% NA

wm_coord_her <- predict_mcod %>%
  group_by(year) %>%
  drop_na(biomass_her) %>% 
  summarise(wlat = weighted.mean(x = lat, w = biomass_her),
            wlon = weighted.mean(x = lon, w = biomass_her)) %>%
  ungroup() %>% 
  mutate(species = "Herring")
#> group_by: one grouping variable (year)
#> drop_na (grouped): removed 21,502 rows (9%), 227,951 rows remaining
#> summarise: now 27 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'species' (character) with one unique value and 0% NA

coord_data <- bind_rows(wm_coord_cod, wm_coord_fle, wm_coord_spr, wm_coord_her)


# Plot
palz <- brewer.pal(n = 5, name = "Dark2")
palz[2] <- "gray50"

mean_oxy_temp %>% 
  filter(variable == "Temperature",
         !year == 1994) %>% 
  ggplot(aes(year, value, color = species, linetype = species, size = species)) + 
  geom_point(size = 2.3, alpha = 0.5) +
  guides(size = "none", linetype = "none") +
  scale_linetype_manual(values = c(1,2,1,1,1)) +
  scale_size_manual(values = c(1,2.3,1,1,1)) +
  scale_color_manual(values = palz, name = "") + 
  stat_smooth(method = "gam", formula = y ~ s(x, k = 3), se = FALSE) +
  labs(x = "Year", y = "Temperature [°C]") + 
  theme(legend.position = "bottom",
        aspect.ratio = 1) +
  NULL
#> filter: removed 86 rows (40%), 130 rows remaining

ggsave("figures/supp/density/weighted_mean_temp_species.png", width = 6.5, height = 8.5, dpi = 600)

# Now do flounder
ggplot(filter(coord_data, !year == 1993), aes(wlon, wlat, color = year)) + 
  facet_wrap(~species, scales = "free") +
  geom_point() +
  geom_path() +
  scale_color_viridis()
#> filter: removed 4 rows (4%), 104 rows remaining

coord_data %>% filter(year %in% c(1994, 2019)) %>% 
  ggplot(aes(wlon, wlat, color = factor(year))) + 
  facet_wrap(~species, scales = "free") +
  geom_point() +
  geom_path() +
  scale_color_viridis(discrete = TRUE)
#> filter: removed 100 rows (93%), 8 rows remaining
#> geom_path: Each group consists of only one observation. Do you need to adjust
#> the group aesthetic?
#> geom_path: Each group consists of only one observation. Do you need to adjust
#> the group aesthetic?
#> geom_path: Each group consists of only one observation. Do you need to adjust
#> the group aesthetic?
#> geom_path: Each group consists of only one observation. Do you need to adjust
#> the group aesthetic?

LongLatToUTM <- function(x, y, zone){
  xy <- data.frame(ID = 1:length(x), X = x, Y = y)
  coordinates(xy) <- c("X", "Y")
  proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")  ## for example
  res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
  return(as.data.frame(res))
}

utm_coordz <- LongLatToUTM(zone = 33, coord_data$wlon, coord_data$wlat)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition

coord_data$wlat_y <- utm_coordz$Y
coord_data$wlon_x <- utm_coordz$X

plot_map_raster +
  facet_wrap(~species) +  
  geom_point(data = filter(coord_data, year %in% c(1994, 2019)),
             aes(wlon_x, wlat_y, color = factor(year)), inherit.aes = FALSE, size = 3) +
  theme_plot() +
  scale_color_brewer(palette = "Dark2", name = "Year") + 
  theme(legend.position = "bottom",
        plot.margin = unit(c(0, 0.1, 0, 0.2), "cm")) + 
  NULL
#> filter: removed 100 rows (93%), 8 rows remaining

plot_map_raster +
  facet_wrap(~species) +  
  geom_point(data = filter(coord_data, !year == c(1993)),
             aes(wlon_x, wlat_y, color = year), inherit.aes = FALSE, size = 1) +
  geom_path(data = filter(coord_data, !year == c(1993)),
             aes(wlon_x, wlat_y, color = year), inherit.aes = FALSE, size = 0.1) +
  theme_plot() +
  scale_color_viridis() + 
  theme(legend.position = "bottom",
        plot.margin = unit(c(0, 0.1, 0, 0.2), "cm"),
        legend.text = element_text(angle = 30)) + 
  xlim(1.8*xmin2, 0.85*xmax2) +
  ylim(1.023*ymin2, 0.976*ymax2) +
  NULL
#> filter: removed 4 rows (4%), 104 rows remaining
#> filter: removed 4 rows (4%), 104 rows remaining
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.
#> Scale for 'y' is already present. Adding another scale for 'y', which will
#> replace the existing scale.

ggsave("figures/supp/density/weighted_mean_coord_species.png", width = 6.5, height = 8.5, dpi = 600)

# pred_grid %>%
#   filter(year > 2010 & sub_div == 25 & depth > 48 & depth < 52) %>% 
#   group_by(year) %>% 
#   summarise(mean_oxy = mean(oxy))
#   
# 
# pred_grid %>%
#   filter(year > 2010 & sub_div == 25 & depth > 50 & depth < 65) %>% 
#   group_by(year) %>% 
#   summarise(mean_oxy = mean(oxy))

# VR combined plot
p1 <- plot_map_raster +
  facet_wrap(~species) +  
  geom_point(data = filter(coord_data, !year == c(1993)),
             aes(wlon_x, wlat_y, color = year), inherit.aes = FALSE, size = 1) +
  geom_path(data = filter(coord_data, !year == c(1993)),
             aes(wlon_x, wlat_y, color = year), inherit.aes = FALSE, size = 0.1) +
  theme_plot() +
  scale_color_viridis() + 
  theme(legend.position = "bottom",
        legend.spacing.x = unit(0, 'cm'),
        legend.spacing.y = unit(0, 'cm'),
        legend.margin=margin(0,0, 0, 0),
        legend.box.margin=margin(-20,-20,-20,-20),
        #strip.placement = NULL,
        #plot.margin = unit(c(0, 0.1, 0, 0.2), "cm"),
        legend.text = element_text(angle = 30)#,
        #plot.margin = unit(c(1, 1, 1, 1), "cm")
        ) + 
  xlim(1.8*xmin2, 0.85*xmax2) +
  ylim(1.023*ymin2, 0.976*ymax2) +
  NULL
#> filter: removed 4 rows (4%), 104 rows remaining
#> filter: removed 4 rows (4%), 104 rows remaining
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.
#> Scale for 'y' is already present. Adding another scale for 'y', which will
#> replace the existing scale.

p2 <- mean_oxy_temp %>% 
  filter(variable == "Temperature",
         !year == 1994) %>% 
  ggplot(aes(year, value, color = species, linetype = species, size = species)) + 
  geom_point(size = 2.3, alpha = 0.5) +
  guides(size = "none", linetype = "none", color = guide_legend(nrow = 1)) +
  scale_linetype_manual(values = c(1,2,1,1,1)) +
  scale_size_manual(values = c(1,2.3,1,1,1)) +
  scale_color_manual(values = palz, name = "") + 
  stat_smooth(method = "gam", formula = y ~ s(x, k = 3), se = FALSE) +
  labs(x = "Year", y = "Temperature [°C]") + 
  theme(legend.position = "bottom",
        legend.text = element_text(size = 6),
        legend.spacing.x = unit(0, 'cm'),
        legend.spacing.y = unit(0, 'cm'),
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(-40,-40,-40,-40),
        aspect.ratio = 1#,
        #plot.margin = unit(c(1, 1, 1, 1), "cm")
        ) +
  NULL
#> filter: removed 86 rows (40%), 130 rows remaining

p1 | p2 

ggsave("figures/supp/density/VR_weighted_mean_coord_species.png", width = 6.5, height = 8.5, dpi = 600)

Marginal effect of oxygen

# Prepare prediction data frame
d2 <- d %>% drop_na(oxy)
nd_oxy <- data.frame(oxy = seq(min(d2$oxy), max(d2$oxy), length.out = 100))

nd_oxy <- nd_oxy %>%
  mutate(year = 2003L,
         depth_sc = 0,
         oxy_sc = (oxy - mean(oxy))/sd(oxy),
         temp_sc = 0)

# Predict
p_margin_oxy <- predict(m_1, newdata = nd_oxy, se_fit = TRUE, re_form = NA)

mar_oxy <- ggplot(p_margin_oxy, aes(oxy, exp(est),
  ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
  geom_line() +
  geom_ribbon(alpha = 0.4) +
  coord_cartesian(ylim = c(30, 220)) +
  theme(aspect.ratio = 1) +
  labs(x = expression(Oxygen[haul]))

Marginal effect of depth

# Prepare prediction data frame
nd_dep <- data.frame(depth = seq(min(d2$depth), max(d2$depth), length.out = 100))

nd_dep <- nd_dep %>%
  mutate(year = 2003L,
         depth_sc = (depth - mean(depth))/sd(depth),
         oxy_sc = 0,
         temp_sc = 0)

# Predict
p_margin_dep <- predict(m_1, newdata = nd_dep, se_fit = TRUE, re_form = NA)

mar_depth <- ggplot(p_margin_dep, aes(depth, exp(est),
  ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
  geom_line() +
  geom_ribbon(alpha = 0.4) +
  coord_cartesian(ylim = c(30, 220)) +
  theme(aspect.ratio = 1) +
  labs(x = expression(Depth[haul]))

Marginal effect of temperature

# Prepare prediction data frame
nd_temp <- data.frame(temp = seq(min(d$temp), max(d$temp), length.out = 100))

nd_temp <- nd_temp %>%
  mutate(year = 2003L,
         depth_sc = 0,
         oxy_sc = 0,
         temp_sc = (temp - mean(temp))/sd(temp))

# Predict
p_margin_temp <- predict(m_1, newdata = nd_temp, se_fit = TRUE, re_form = NA)

mar_temp <- ggplot(p_margin_temp, aes(temp, exp(est),
  ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
  geom_line() +
  geom_ribbon(alpha = 0.4) +
  coord_cartesian(ylim = c(30, 220)) +
  theme(aspect.ratio = 1) +
  labs(x = expression(Temp[haul]))

Plot together

# mar_temp + mar_depth + mar_oxy +
#   plot_annotation(tag_levels = 'A') 
# 
# ggsave("figures/supp/density/marginal_effects.png", width = 6.5, height = 6.5, dpi = 600)

p_margin_dep2 <- p_margin_dep %>%
  mutate(variable = "Depth") %>% 
  dplyr::select(est, est_se, depth_sc, variable) %>% 
  rename(value = depth_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_oxy2 <- p_margin_oxy %>%
  mutate(variable = "Oxygen") %>% 
  dplyr::select(est, est_se, oxy_sc, variable) %>% 
  rename(value = oxy_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

p_margin_tem2 <- p_margin_temp %>%
  mutate(variable = "Temperature") %>% 
  dplyr::select(est, est_se, temp_sc, variable) %>% 
  rename(value = temp_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)

margs <- bind_rows(p_margin_dep2, p_margin_oxy2, p_margin_tem2)

ggplot(margs, aes(value, exp(est), ymin = exp(est - 1.96 * est_se),
                  ymax = exp(est + 1.96 * est_se))) +
  geom_ribbon(alpha = 0.4) +
  geom_line() +
  facet_wrap(~variable, scales = "free_x") +
  labs(y = "Biomass density",
       x = "Scaled value") +
  theme_plot() +
  theme(axis.text.x = element_text(angle = 0)) +
  NULL


ggsave("figures/supp/density/marginal_effects.png", width = 6.5, height = 6.5, dpi = 600)
knitr::knit_exit()